library("tidyverse")
library("firatheme")
library("estimatr")

# load data
load("~/df_replication.RData")

model2 <- lm_robust(click ~ issue, data = df_replication, se_type="HC2")
model4 <- lm_robust(click ~ issue + region + party_name, data = df_replication, se_type="HC2")
model6 <- lm_robust(click ~ cue*issue + region + party_name, data = df_replication, se_type="HC2")


figure9.1 <- tibble(
  Treatment = c("Climate Change","Public Schools"),
  avg = c(round(coef(model2)[1]*100,1), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1)), 
  lower = c(round(coef(model2)[1]*100,1) - round(model2[["std.error"]][["(Intercept)"]]*100,4), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1) - round(model2[["std.error"]][["issueschool"]]*100,4)),
  upper = c(round(coef(model2)[1]*100,1) + round(model2[["std.error"]][["(Intercept)"]]*100,4), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1) + round(model2[["std.error"]][["issueschool"]]*100,4)))

ggplot(figure9.1, aes(x= Treatment,y=avg, fill=Treatment)) +
  geom_bar(stat="identity",position="dodge", color= "black", width=.25)+
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(.7))+
  geom_text(aes(y=upper, label=round(avg,3), size = 12), position=position_dodge(width=0.7), vjust=-.75)+
  labs(y="Click Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 35)) + theme_fira() +  ggtitle("a") +
  theme(text = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = 0))


estimates <- c(0,round(coef(model4)[2]*100,1))
lower95 <- c(0,confint(model4)[2,1] * 100) 
upper95 <- c(0,confint(model4)[2,2] * 100) 
lower90 <- c(0,confint(model4, level = 0.90)[2,1] * 100) 
upper90 <- c(0,confint(model4, level = 0.90)[2,2] * 100) 
labels <- c("climate change", "public \nschools")
plot9.2 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot - problem
ggplot(plot9.2, aes(x=labels, y=estimates, color=labels)) + #color = type 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) +
  scale_color_manual(values=c("gray50", "gray50")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "Note: 'Climate Change' is reference category. Error bars represent 90% and 95% Confidence Intervals") +
  ggtitle(c("b")) +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = 0), strip.text.x = element_blank())


a <- prediction(model6, at = list(issue = "school", cue = "problem_indicators"))
b <- prediction(model6, at = list(issue = "climate", cue = "problem_indicators"))
c <- prediction(model6, at = list(issue = "school", cue = "public_opinion"))
d <- prediction(model6, at = list(issue = "climate", cue = "public_opinion"))
e <- prediction(model6, at = list(issue = "school", cue = "rival_parties"))
f <- prediction(model6, at = list(issue = "climate", cue = "rival_parties"))
g <- prediction(model6, at = list(issue = "school", cue = "the_media"))
h <- prediction(model6, at = list(issue = "climate", cue = "the_media"))


figure9.3 <- tibble(
  Treatment = c("Problem Indicators","Problem Indicators","Public Opinion", "Public Opinion", "Party Positions", "Party Positions", "Media Reporting", "Media Reporting"),
  Problem = c("Public Schools","Climate Change","Public Schools", "Climate Change", "Public Schools", "Climate Change", "Public Schools", "Climate Change"),
  avg = c(mean(a$fitted)*100, mean(b$fitted)*100, mean(c$fitted)*100, mean(d$fitted)*100, mean(e$fitted)*100, mean(f$fitted)*100, mean(g$fitted)*100, mean(h$fitted)*100), 
  lower = c(mean(a$fitted)*100 - mean(a$se.fitted)*100, mean(b$fitted)*100 - mean(b$se.fitted)*100, mean(c$fitted)*100 - mean(c$se.fitted)*100, mean(d$fitted)*100 - mean(d$se.fitted)*100, mean(e$fitted)*100 - mean(e$se.fitted)*100, mean(f$fitted)*100 - mean(f$se.fitted)*100, mean(g$fitted)*100 - mean(g$se.fitted)*100, mean(h$fitted)*100 - mean(h$se.fitted)*100),
  upper = c(mean(a$fitted)*100 + mean(a$se.fitted)*100, mean(b$fitted)*100 + mean(b$se.fitted)*100, mean(c$fitted)*100 + mean(c$se.fitted)*100, mean(d$fitted)*100 + mean(d$se.fitted)*100, mean(e$fitted)*100 + mean(e$se.fitted)*100, mean(f$fitted)*100 + mean(f$se.fitted)*100, mean(g$fitted)*100 + mean(g$se.fitted)*100, mean(h$fitted)*100 + mean(h$se.fitted)*100))

ggplot(figure9.3, aes(x= reorder(Treatment,-avg), y=avg, fill=Problem)) +
  geom_bar(stat="identity", position="dodge", color= "black", width=.25) +
  geom_linerange(aes(ymin=lower, ymax=upper), position= position_dodge(.25)) +
  geom_text(aes(y=upper, label=round(avg,1)), size = 3, position=position_dodge(width=0.7), vjust=-1.75)+
  labs(y="Click Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50", "grey80","gray50", "grey30","gray50", "grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 40)) + theme_fira() + ggtitle("c") +
  theme(text = element_text(size=8), legend.position= "bottom", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = 0))

